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Abstract 

Visually induced neuronal activity in VI displays a marked gamma-band component which is 
modulated by stimulus properties. It has been argued that synchronized oscillations contribute 
to these gamma-band activity. However, analysis of Local Field Potentials (LFPs) across dif- 
ferent experiments reveals considerable diversity in the degree of oscillatory behavior of this 
induced activity. Contrast-dependent power enhancements can indeed occur over a broad band 
in the gamma frequency range and spectral peaks may not arise at all. Furthermore, even when 
oscillations are observed, they undergo temporal decorrelation over very few cycles. This is not 
easily accounted for in previous network modeling of gamma oscillations. 

We argue here that interactions between cortical layers can be responsible for this fast decor- 
relation. Wc study a model of a VI hypercolumn, embedding a simplified description of the 
multi-layered structure of the cortex. When the stimulus contrast is low, the induced activity 
is only weakly synchronous and the network resonates transiently without developing collec- 
tive oscillations. When the contrast is high, on the other hand, the induced activity undergoes 
synchronous oscillations with an irregular spatiotemporal structure expressing a synchronous 
chaotic state. As a consequence the population activity undergoes fast temporal decorrela- 
tion, with concomitant rapid damping of the oscillations in LFPs autocorrelograms and peak 
broadening in LFPs power spectra. 

We show that the strength of the inter-layer coupling crucially affects this spatiotemporal 
structure. We predict that layer VI inactivation should induce global changes in the spectral 
properties of induced LFPs, reflecting their slower temporal decorrelation in the absence of inter- 
layer feedback. Finally, we argue that the mechanism underlying the emergence of synchronous 
chaos in our model is in fact very general. It stems from the fact that gamma oscillations induced 
by local delayed inhibition tend to develop chaos when coupled by sufficiently strong excitation. 

Author Summary 

Visual stimulation elicits neuronal responses in visual cortex. When the contrast of the used 
stimuli increases, the power of this induced activity is boosted over a broad frequency range (30- 
100 Hz), called the "gamma band". It would be tempting to hypothesize that this phenomenon 
is due to the emergence of oscillations in which many neurons fire collectively in a rhythmic 
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way. However, previous models trying to explain contrast-related power enhancements using 
synchronous oscillations failed to reproduce the observed spectra, because they originated un- 
realistically sharp spectral peaks. The aim of our study is to reconcile synchronous oscillations 
with broad-band power spectra. We argue here that, thanks to the interaction between neu- 
ronal populations at different depths in the cortical tissue, the induced oscillatory responses 
are synchronous, but, at the same time, chaotic. The chaotic nature of the dynamics makes 
possible to have broad-band power spectra together with synchrony. Our modeling study allows 
us formulating qualitative experimental predictions, that provide a potential test for our theory. 
We predict that if the interactions between cortical layers are suppressed, for instance by inac- 
tivating neurons in deep layers, the induced responses might become more regular and narrow 
isolated peaks might develop in their power spectra. 



Introduction 



An increase of activity in the gamma band (30-100 Hz) is observed in Local Field Potential (LFP) 
and Multi-Unit Activity (MUA) recordings [T|-|14], as well as in EEG and Electrocorticogram 



studies [15,16 in primary visual cortex (VI) upon visual stimulation. Gamma activity is mod' 
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ulated by properties of the presented stimulus, such as orientation [2, 14, 17 
velocity]^ 
bands [19|[20 



or size [12] , much more strongly than the change in power in other frequency 
. Local GABA-ergic interneuronal networks are thought to play a key role in the 
production of neuronal activity in the gamma range ( [21], see ^22j for a review), as upheld as 



well by recent results obtained through optogenetic techniques in- vivo [23 , 24 



Modeling works have provided a theoretical basis to account for the way in which networks 
of inhibitory interneurons can generate synchronous oscillatory activity in the gamma range 
[25| - |30[ . In brief, in one possible scenario, the dynamics of the inhibitory post-synaptic currents 
is non-instantaneous (due to axonal delays, but also simply to finite synaptic time-constants). 
This contributes to create narrow time-windows in which excitatory and inhibitory neurons 
can fire closely in-phase, before being prevented to do so by a delayed inhibitory feedback. 
Therefore delayed inhibition, without need of an active involvement of excitatory populations, 
is capable inducing collective synchronous oscillations in neuronal activity. The frequency of 
these oscillations falls in the gamma band if the synaptic time constant of the inhibition is in an 
appropriate range. If a network operates in such a synchronous regime the neurons are engaged 
into approximately periodic collective oscillations involving a macroscopically large number of 
neurons. Therefore these oscillations are weakly affected by local noise and they maintain 
coherence over arbitrarily long time intervals. Power spectra of population observables of the 
network activity (e.g. LFP or MUA) exhibit narrow harmonic-like peaks and the damping of 
the corresponding autocorrelograms is slow. 

Peaks in the gamma-band have been identified in the LFP or MUA spectra of induced 
activity in- vivo in VI [l}|4, 12 . However, in general these peaks are very broad and in many 
cases they are virtually indistinguishable as the stimulus-modulated gamma power of the signals 
spreads across a broad-band frequency interval [7[[9| [Io[[3l}|33] . Characterization of the spatio- 
temporal structure of the gamma induced activity by means of auto-correlations (AC) and 
cross-correlations (CC) of single-unit, multi-unit and LFP signals has also revealed that the 
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neuronal activity has a tendency to oscillate, which can be stronger or weaker, depending on the 
considered experiment. In some cases the oscillatory components of AGs and CCs of the induced 
activity display many cycles before getting damped [l|[2, 12, 14 . In other cases, however, the 



oscillations are completely damped after one or two cycles [3|[8 ,13 ,17 . The existence of different 
dynamical regimes might underlie this observed diversity. 

For the mathematical abstraction of infinitely large networks, sharp boundaries between 
asynchronous and synchronous dynamical states exist 134] , but for networks of a finite size such 



transitions are fuzzier 25 , 34 , 35 . Consequently, if the network does not operate too far from the 
instability to collective oscillations, in a regime which is formally defined as asynchronous — see 
[34p5 and below for the definition — , the dominant normal modes of the network, which describe 



its response to small perturbations, can display damped oscillations at gamma frequencies. 
Local noise can excite these modes, inducing short-lived episodes of synchronous oscillatory 
activity. However, since these episodes are transient, the subsequent increase in power at gamma 
frequencies is broad-band, induced broad band gamma power increases in VI can therefore be 
accounted for if one assumes that the VI network operates in such an asynchronous regime at the 
edge of developing synchrony 36-38 . In this regime, correlations in the spikes as well as in the 
membrane potentials of pairs of neurons are in general weak unless the neurons are connected via 
strong and direct synapses. However, in order to get a significant, although damped, oscillatory 
component in the macroscopic activity, the network must be "at the edge of synchronization". 
Parameters have to be tuned in such a way to be close to an instability toward fully- developed 
synchronous oscillations, and this tuning have to be tighter, the larger the size of the recruited 



network 25 , 29 . It is not clear how the required fine tuning would be satisfied given the range 
of experimental conditions in which gamma oscillations have been observed. 

In the present study we explore another scenario which reconciles collective synchronous 
activity with broad-band spectral modulations and robust fast decoherence. It is based on a 
mechanism proposed recently for the emergence of synchronous chaos in recurrent neural net- 



works [39-41 . In this mechanism, clusters of neuron undergo a synchronous gamma oscillation 
due to local mutual inhibition. These collective gamma oscillations become chaotic when the 
neuronal clusters are allowed to interact through longer-range excitation. The resulting overall 
patterns of activity are characterized by synchrony at the population level, but at a same time 
display a characteristic lack of temporal regularity due to chaos. As a consequence, the power of 
this activity spreads over a broad interval of frequencies and the oscillatory components of the 
autocorrelograms of neuronal activity and LFP signals are rapidly damped within a few tenths 
of a millisecond. In this alternative regime, correlations in the spikes of pairs of neurons are 
still weak and go together with the sparseness of the firing, but correlations in their membrane 
potentials can be strong. 

We present here a model of a hypercolumn in VI, endowed with a simplified multi- layer 
architecture. In order to explain broad-band contrast-dependent spectral modulations in terms 
of synchronous chaos, we need to identify distinct interacting oscillators within the local cortical 
circuit. We hypothesize that neuronal populations within different thalamo-recipient cortical 
layers are set into oscillation by increased driving and that the mutual interaction between 
these populations, mediated by inter-layer synaptic connections, supports the development of 
synchronous chaos. This hypothesis is backed up by anatomical evidence. Thalamo-cortical 
synapses, providing direct sensory-induced driving, indeed target cortical layer IV but also, to a 
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lesser extent, layer VI 142-46 



within each thalamo- recipient layer 145 47 , 48 



Extensive networks of recurrent inhibitory connections are present 
supporting local generation of oscillations at 



multiple depths in the cortical tissue. Finally, stereotyped circuit motifs provide a bidirectional 



poly-synaptic connection loop between thalamo-recipient layers 44-46,49-51 



Relying on extensive numerical simulations, we show that our model displays broad-band 
gamma modulations of the spectra of LFPs upon stimulation of the network at low as well 
as at high contrast. Whereas this induced activity is asynchronous at low contrast, it develops 
synchrony on a macroscopic scale when the contrast increases. Therefore we argue that the broad 
band gamma power observed in recorded LFP spectra in VI is compatible with the existence of 
visually induced synchronous oscillatory neuronal dynamics. 



Results 

Multi-layer hypercolumn model 

We model a functional hypercolumn in primary visual cortex as a large recurrent network of 
spiking integrate-and-fire-type neurons. To account in a simplified way for the layered structure 
of the visual cortex — a cartoon of which is shown in Figure [T]A — the model network consists of 
two sub-networks, schematically representing layers I to IV and layers V to VI. We denote these 
two sub- networks as the upper and lower layer respectively (Figure [Tj3). Each of these layers 
comprises Ne excitatory and Ni inhibitory neurons, for a total number of = 2{Ne + Nj) 
neurons in the network. Most of the simulations in this study are performed taking Ne = 4000 
excitatory and Ni = 1000 inhibitory neurons per layer, leading a total of = 10000 neurons 
in the model hypercolumn. This number is one order of magnitude smaller than estimates of 
the number of neurons in a real VI hypercolumn based on neuronal densities recently measured 



by 52 . However, it leads to dynamical behaviors similar to larger network sizes (see following 
scaling analyses) and constitutes a compromise for efficient and fast simulations. 

Each layer is described by a network with the geometry of a ring as depicted in Figure [Tp, 



with neurons labeled by angular coordinates, i}, ranging from -90 to -|-90 degrees 53 ,54|. The 
connections between neurons within each layer are random, with connection probabilities that 
depend on the angular distance between pre- and post-synaptic neurons. Spatial averages and 
spatial modulations of connection probabilities are set independently for the various kinds of 
connections (e.g. excitatory-to-excitatory, excitatory-to-inhibitory, inhibitory-to-excitatory or 
inhibitory-to- inhibitory) , thus making it possible to vary the spatial profiles of net synaptic in- 
teractions (see Figure[lp, E, F). Excitatory and inhibitory inter-layer connections are also ran- 
dom and spatially modulated. All the external inputs to the network are modeled as stochastic 
processes (see Methods section). The neurons receive an external non-selective noisy current 
representing background inputs to VI from other cortical and subcortical areas and a weakly 
tuned noisy current which represents visually induced inputs to VI from converging Lateral 
Geniculate Nucleus (LGN) synapses [55j. Note that the two main thalamo-recipient layers, i.e. 
layers VI and IV, are embedded within two distinct model layers. 

Our two- layer circuit embeds in a simplified manner several known features of the stereotyp- 
ical interlaminar anatomy of the columnar microcircuit, in particular, the existence of a layers 



IV to VI to IV feedback loop 44 ,46 ,50 . Furthermore, a different degree of spatial modulation 
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for inter-layer excitation and inhibition mimic the on-center off-surround arrangement of layers 
VI to IV projections [56]. In the simulations described below we assume that the LGN input to 
the lower layer is weaker (by a factor of 2) than the input to upper layer to account for the fact 
that thalamo-cortical synapses reaching layer VI are smaller in number than those reaching layer 



IV 45 . We also assume that latencies for inter-layer connections are longer than for intra-layer 
connections, thus accounting for the multisynaptic nature of this coupling. Our assumptions on 
the connectivity, external inputs and latencies are further commented upon in the Discussion 
section. 

In order to analyze the role of the interlayer interactions in shaping the spatiotemporal 
dynamics of our model hypercolumn, we introduce a parameter < F < 1 which homogeneously 
rescales the strength of excitatory and inhibitory connections between layers. For F = 1 the 
interactions between the layers assume their maximum strength. For F = the layers are 
completely independent. In the following, we consider first the dynamics of the network at full 
coupling strength, F = 1. 

Orientation tuning and contrast dependence of induced response 

In absence of "visual" stimuli (contrast level C = 0%), the model hypercolumn is driven only by 
the non-selective background input. The resulting spontaneous activity is heterogeneous across 
the neurons with average firing rates of 1.2 it 0.4 Hz and 5 it 3 Hz for excitatory and inhibitory 
neurons, respectively. Differences in the spontaneous firing rate distributions for upper and 
lower layers are not statistically significant at the 5% confidence level. The spontaneous firing of 
the neurons is highly irregular due to the stochasticity of the inputs. For instance, the average 
coefficient of variation (CV) of the interspike histogram of excitatory and inhibitory neurons in 
the upper layer is CV = 0.9 it 0.2. More details about rate and CV distributions can be found 
in Appendix SI. 

The profile of the activity induced by an oriented stimulus in both layers, is localized and 
centered at an angular coordinate corresponding to the stimulus orientation. Hence, the neuronal 
responses are selective to the stimulus orientation. The tuning curves of individual neurons 
display some heterogeneity in their broadness, as exhibited by distributions of peak response 
rates, circular variance and skewness of the tuning curves (reported in Appendix S3). 

Figure [2]A displays the population average tuning curve for various contrast levels for ex- 
citatory neurons in the upper layer. Comparison between tuning curves at different contrasts 
reveals that tuning width is approximately contrast invariant and that the larger deviations are 
observed for small contrast levels (tuning curves normalized to the peak are plotted in Appendix 



S4) . This invariance is achieved as an effect of noise in synaptic inputs [57 , 58 . 

The preferred responses of the excitatory neurons vary non-linearly with the contrast as 
depicted in Figure [2^, where the population average Contrast Response Functions (CRFs) are 
plotted for excitatory neurons in the upper layer. It can be fitted by an hyperbolic ratio function 
(see Methods section), with mid-range contrast C50 ~ 15% and an exponent of n 5 (upper 
layer neurons). This nonlinear dependence stems from the fact that increased sensory-driving 
yields larger inhibitory neurons activity which in turn is responsible for the saturation of the 
excitatory population response ^59j. The CRFs of inhibitory neurons show a much weaker 
tendency to saturation at large stimulus contrasts which is due to the logarithmic dependency 
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on the contrast of their external input. The CRFs of single neurons are heterogeneous, in 
qualitative agreement with experimental reports [60| (see Appendix S3). The contrast response 
functions of the lower layer are homologous, but the induced responses are approximately twofold 
smaller, due to the weaker LGN driving. 



The dynamical state of the network depends on the stimulus contrast 

For zero contrast, the synchrony level in the spontaneous neuronal activity is small, as denoted 
by a small value of the synchrony factor x- This factor, defined in the Methods section, quantifies 
global synchrony over a network and is bounded between and 1. For a network of size N = 
10000, the synchrony factor for spontaneous activity assumes the value x = 0.02. Furthermore, 
as shown in Figure 6B (grey line), it vanishes consistently as x oc for larger network sizes, 

allowing us to classify formally the (asymptotic) state of the network as "asynchronous" (see 
Methods section and |34[|35| ). 

The single neuron and population responses of the network to a low stimulus contrast, C = 
2% is illustrated in Figure [3} The raster plot of the spike activity of all the excitatory neurons in 
the upper layer is plotted in Figure|3|\. It suggests that the firing is highly irregular (the mean CV 
of the upper layer excitatory neurons is 0.9ib0.1, see Appendix SI) and that the network activity 
of the network is only weakly synchronized. This is confirmed in Figure [3^ where the spike trains 
of six upper layer cells stimulated within ±5° from their preferred orientation are plotted. The 
neurons fire without any noticeable synchrony. Figure [3p displays the voltage traces of two of 
these neurons. The comparison between the sub-threshold fiuctuations in the two traces does 
not reveal any significant correlation. To further quantify the correlations in the supra and 
subthreshold activity of the neurons we compute the zero delay pairwise correlation coefficients 
(CCos) of the spikes and the membrane potential traces for a large number of pairs formed by 
highly active neurons with preferred orientation within ±9° from the presented stimulus (see 
Methods section and Appendix S2 for details). The resulting histograms are shown in Figure [3|D 
(spikes: left, cyan color; voltage: right, blue color). They are peaked around zero with a mean 
statistically indistinguishable from zero (0.000 it 0.001 for spikes and voltage). Almost all the 
CCos are weak for the spikes as well as for the voltage traces (CCos larger than 0.25 occur only 
for 2% of the pairs when considering spike CCos, and for 0.1% of the pairs when considering 
voltage trace CCos). These results are consistent with a very weak synchronization in the 
network activity. This is in line with the small value of the synchronization factor, which is only 
X = 0.03. Auto- and crosscorrelograms of spike trains and membrane potential traces of three 
representative neurons are also shown in Figure [5)^,B. The pairwise crosscorrelograms of both 
spikes and voltages do not display any persistent oscillatory component, even when two cells 
share a same orientation preference. 

The dynamical state of the network is qualitatively different for a high contrast stimulus. For 
C = 95% the neurons are engaged into a collective pattern of synchronous oscillations in contrast 
to what happens for C = 2%. This is clear from the raster plot in Figure |4]A. Figure |4p plots the 
membrane potential traces of two neurons. Comparison of these traces suggests that now the 
subthreshold membrane fiuctuations of the neurons are strongly correlated across the network. 
As a matter of fact, the synchrony factor, Xi which characterizes the degree of synchrony in 
the subthreshold activity at the network level, is x = 0.75. However action potentials are 
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much less synchronized, as suggested by the comparison of the spike trains of the six neurons 
plotted in Figure [4^: although multi- neuron coincidences in firing (denoted by vertical grey 
bars) can be detected, the overall synchrony is weak. This substantial difference in the strength 
of the pair correlations in supra and subthreshold activities is clear in Figure |4p. All the CCos 
of the subthreshold membrane potentials (red histogram) are large and sharply distributed 
around 0.75 (standard deviation of ±0.03) whereas the distribution of the spike trains CCos 
(magenta histogram) has a mean which is only 0.097 it 0.002. Remarkably, the firing activity 
continues to be highly irregular, despite the high degree of synchrony (mean CV of upper layer 
excitatory neurons is CV= 1.0 it 0.2, see Appendix SI). Auto- and crosscorrelograms of spike 
trains and membrane potential traces of three representative neurons are shown in Figure [5p,D. 
The pairwise crosscorrelograms of voltages display now a clear oscillatory structure, which is 
however completely damped after only two or three cycles. Note that oscillatory correlations 
are evident even when the difference of preferred orientation is large (> 20°). Note as well that 
pairwise crosscorrelograms of spike trains do not display any marked oscillation even when the 
two considered cells have similar preferred orientations. We stress that the small mean value 
CCos and the lack of a clear oscillatory structure in the crosscorrelograms for spike trains, in 
both the low and the strong contrast case, is associated to the irregularity and the sparseness 
of single neuron firing. 

These results indicate that synchrony in the population activity increases with the contrast. 
As a matter of fact, the synchrony measure x varies abruptly around a contrast value of ~10%, 
as shown in Figure [6]A. This is even sharper with larger network sizes (compare in Figure [gJB, 
the solid line which is for = 10000 with the dashed line which is for N = 40000). Moreover, a 
systematic analysis of the dependency of x o^i the size N reveals that for C < 10%, x {low con- 
trast re5'ime)vanishes consistently with N, x while for C > 10% {large contrast regime) it 
converges toward a constant non zero value (Figure [6^). Hence, the network operates in quali- 
tatively different regimes at low and high contrast. Whereas the network state can be classified 
as asynchronous in the low contrast regime, it is synchronous in the high contrast regime. This 
sharp variation of synchrony is indicative of a phase-transition occurring for increasing contrast, 
due to an increased drive to the network (see Discussion). 

To characterize further how the population dynamics depend on the contrast we compute 
the autocorrelation, AC{t) of the LFP signals induced by stimuli oriented at the preferred 
orientation of the recording site (see Methods for the way we relate the LFP signals to the 
neuronal activities in the framework of our model and Appendix S6 for examples of LFP traces) . 
The result for low contrast, C = 2%, is plotted in Figure [T]^, B. The amplitude of the (non- 
normalized) AC at zero delay, AC{0), is small and decreases with the network size as 1/^/N. 
Similarly, the small oscillatory component of the AC disappears gradually for increasing network 
sizes (Figure [T^). This is because the network state is asynchronous and in a larger network 
more cells contribute to the LFP signal (see also Methods section). 

The fact that at high contrast, C = 95%, the network is engaged in collective synchronous 
activity is manifest in Figure [7p,D: AC{0) is now large and it does not vanish in the large limit 
and is almost independent of A^ for > 10000. However, and remarkably, the induced dynamics 
exhibit a spatio-temporal structure which is more complex than a periodic regular oscillation of 
the population activity: the time interval between consecutive episodes of synchronous activity 
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displays cycle-to-cycle fluctuations as can be observed in the raster plotted in Figure |4|A). As a 
result, the LFP autocorrelogram is rapidly damped. Although it displays some secondary peaks 
their amplitudes are very small as shown in Figure [7p. The damping of the AC oscillations is even 
faster for larger network sizes (Figure [7p). Note that autocorrelations for intermediate contrast 
values are also rapidly damped (see Appendix S6). A moderate tendency to period doubling, 
manifested by a second autocorrelogram peak slightly larger than the first autocorrelogram peak, 
has not been reported experimentally. We remark however that this is an accidental feature, 
which is no more observed, for instance, for larger network sizes or stronger inter-layer coupling. 

LFPs induced by non-preferred stimulus directions display as well oscillatory components, for 
both low and high contrasts. Induced LFPs are correlated over the entire ring network as revealed 
by crosscorrelation analysis, confirming that sub-threshold coherence can exist independently 
from correlations in spiking activity (see Appendix S6). 

Finally, we consider the spectral properties of induced LFPs, and their relation with MUA 
observed at a same location. The dependency on the contrast of the power spectra of the 
LFPs induced by preferred-orientation stimuli is shown in Figure |8j^. The low-frequency part 
of the power spectra is weakly dependent on the stimulus contrast. Rather, it is shaped by the 
properties of cortical background activity, modeled as a stochastic Ornstein-Uhlenbeck noise with 
a frequency cutoff (see Methods section and [6l|). This should be compared to the boosting of the 
power as the contrast increases for frequencies > 30 Hz. Although the network activity becomes 
much more synchronous at large contrast as explained above, power spectrum modulations 
are not limited to narrow peaks, but, even at the highest contrast, the whole frequency range 
comprised between 30 and 100 Hz is boosted. In this same broad frequency range in which 
contrast-dependent power modulations occur, the LFP displays phase-synchronization with the 
MUA at a same location, as measured by a MUA-LFP coherence increasing with contrast (see 
FigurejS^). Interestingly, the MUA-LFP coherence, even at full contrast, rises only at an average 
peak level of approximately 0.3, compatible with physiologic ranges of synchronization [9,62 



This can be explained by the random-like variability of single neuron firing — inherited by 
the MUA signal, which refiects the spiking activity of only a limited number of single units 
(see Methods section) — , but also by the lack of phase autocoherence in the LFP signal itself 
(cfr. [63]). 

The spatio-temporal structure of the induced activities in the lower and in the upper layers 
are similar. In our simulations, the lower layer average firing rate is approximately half of that 
in the upper layer, reflecting weaker driving from LGN. Cross-correlation analysis of the LFPs 
in the two layers shows that the lower layer oscillations lag behind those in the upper layers 
(see Appendix S5). Note that larger response latencies in deep layers have been experimentally 
observed in specific conditions [64||65| . However, the multi-layer structure in our model is 
too schematic to capture quantitatively such inter-layer relations. In particular, the difference 
in response rate and the exact locking pattern between layers are not robust in our model 
against changes in the parameters of LGN input and inter-layer coupling. On the contrary, 
the synchronization and the fast decorrelation of induced oscillations are robust qualitative 
properties (see later Discussion). 



9 



The role of inter-layer coupling in destroying the temporal coherence of the 
oscillations 

In order to explore the role played by the inter-layer interactions, we investigate in the following 
how the dynamics in the high contrast regime is affected by a change of this coupling. More 
specifically, we rescale the peak conductances of all the synapses between cells in different layers 
by a same factor < F < 1 (here F = 1 and F = correspond respectively to fully coupled and 
fully decoupled layers). 

Upon layer-decoupling the mean firing rate of the excitatory and inhibitory cells increases in 
the upper layer (Figure|9j^). However response rate changes are highly heterogeneous across cells 
and, in some cases, the peak rate is even slightly reduced. An analogous heterogeneity is observed 
in the changes of preferred orientation, skewness and tuning width. However, even though 
changes after complete layer decoupling can be significant for specific cells, the distribution of 
tuning curve parameters over the entire upper layer excitatory neurons population is only weakly 
altered. Details are shown in Appendix S7. 

Another effect of layer decoupling, albeit moderate, is that the degree of synchrony in induced 
activity decreases monotonically with F (Figure [9|3). For instance, the synchrony factor is x = 
0.75 for F = 1, but decreases to x = 0.71 when F = 0.5, and drops further to x = 0-67 for fully 
decoupled layers. 

The most striking consequence of the reduction in inter-layer coupling is the progressive 
qualitative change in the shape of the LFP autocorrelograms and power spectra as F decreases. 
This is depicted in Figure [9j For 80% coupling strength (F = 0.8), the autocorrelogram of LFP 
and the corresponding power spectrum are similar to what is found in the fully-coupled case (fast 
temporal decorrelation and broad plateau- like peak in the gamma spectral band, see Figure |9p, 
D). However, for a 60% coupling strength (F = 0.6), the LFP temporal decorrelation becomes 
considerably slower and the envelope of the autocorrelogram displays amplitude modulations 
indicating that the LFP signal is quasi-periodic. In parallel, the gamma-band spectral plateau 
is replaced by a system of narrow peaks at incommensurate frequencies. The raster plot of 
activity (not shown) continues to display a temporally irregular oscillation; however spatial 
fluctuations in the width of consecutive bumps of spiking activity are reduced with respect to 
the fully-coupled case. For further reduction of the interlayer coupling to F = 0.2, the LFP 
autocorrelogram starts revealing periodicity of the signal over long time scales. The multiple 
narrow spectral resonances merge into a single prominent resonance in the gamma-band and 
secondary harmonic peaks also appear. Finally, for F = (Figure |9^, F), the LFPs are still 
substantially autocorrelated after several hundredths of ms. Spectra in the synchronous regime 
are harmonic at any contrast level. More details about the high contrast regime for completely 
uncoupled layers are presented in Appendix S7. 

Interestingly, qualitative modifications of the population dynamics when F is varied do not 
occur in the low contrast regime, in which collective oscillations do not develop. As a matter of 
fact, independently of the coupling strength F, induced activity is asynchronous. Spiking and 
LFP responses to a low contrast stimulus between completely uncoupled or fully coupled layers 
are practically indistinguishable (not shown). 
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Stimulus repetition and chaotic sensitivity to initial conditions 

Up to now we have focused on the response of the network to a time independent stimulus. Here 
we show that the inter-layer coupling also strongly affects the response of the model hypercolumn 
induced by an external input which varies periodically, representing visual stimuli to VI in the 
form of flashed or drifting gratings. In this situation, we characterize the neuronal responses by 
means of peristimulus time histograms (PSTHs) which express the probability of observing the 
firing of a spike at a given time relative to the onset of each stimulus presentation (see Methods 
section). In the following, we focus on high contrast stimuli. 



The PSTH for F = 1 is shown in Figure 10 A.. At the onset of the stimulus the probability of 
firing increases sharply, followed by a transient phase of reduced firing. This feature is not evident 
in experimental PSTHs. It is due to the strongly synchronous recruitment of recurrent inhibition 
which follows the initial burst of activity, triggered by the rise of external inputs (instantaneous 
in our model). Notwithstanding, after a few tenths of a ms the firing probability rises again 
and remains then almost constant. This reflects the fact that the population responses are 
highly variable across trials as is clear in Figure [TO^. In each trial the response of the network 
consists of a sequence of episodes in which the neurons tend to fire together. However, there 
are substantial trial-to-trial fluctuations in the timing of these episodes and their amplitude (i.e. 
the numbers of recruited cells). Consequently, although the presentations of the stimulus do 
give rise to synchronous activity, the PSTH histogram averaged over many trials is almost flat 
after a peri-stimulus time on the order of the short temporal decorrelation time of the induced 
oscillation. 

In contrast, for fully decoupled layers (F = 0), the PSTH averaged over many trials exhibits a 
long- lasting, although damped population oscillation, as plotted in Figure [TOp. This is because 
when the layers are decoupled the oscillations generated inside the layers are close to being 
periodic and they maintain coherence over several hundred milliseconds. Hence the timing of 
the oscillations does not fluctuate much across trials (Figure [TOp). Population oscillations are 
thus masked by averaging across multiple stimulus repetitions only after many cycles. 

The large trial-to-trial variability displayed by the network for F = 1 (Figure [lOp) indicates 
a strong sensitivity to initial conditions (i.e. the network configuration at the onset of the 
stimulus). To further illustrate this sensitivity, we perturb the dynamics of the system by 
omitting artificially a single spike in a single neuron (out of = 10000) at the center of the bump 
of induced activity and we compare then the perturbed and the unperturbed dynamics. The 
results of this numerical simulation are illustrated by Figure [Tl] As visible from the raster plot 
(Figure [m^) and the population rate histogram (Figure[Tl^) of the upper layer induced activity 
(at full contrast), the perturbed and the unperturbed collective oscillations can be distinguished 
already after one oscillation cycle. After a few cycles, they have completely diverged. Such 
extreme sensitivity to perturbations or initial conditions is strongly indicative of dynamical 
chaos [66] . The sequence of states observed in our model for decreasing F (from irregular to 
quasi-periodic to periodic, see Figure [9]C,D) also suggests that chaos might emerge for strong 
inter-layer coupling and that its onset might occur according to a quasi-periodic scenario |66|67| . 
This is indeed one of the possible scenarios for the transition to chaos occurring in a related rate 
model |41j|. As we discuss in detail in the Appendix SIO, the chaotic nature of the dynamics of 
the network for F = 1 and high contrast stimuli can be assessed by an estimation of its largest 
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Lyapunov exponent Xmax 66 . A positive value of this Lyapunov exponent is the manifestation 



of deterministic chaos, denoting exponentially fast separation of trajectories. Using techniques 
of non-linear time-series analysis 68 applied to very long stationary time-series of LFP from 
our model (see Methods section and Appendix SIO), we obtain the estimate Xmax = 2.2 it 0.6 
ms~^, which is indeed positive. Interestingly, the dynamics of the network with uncoupled layers 
(r = 0) fails to display a positive Lyapunov exponent (see Appendix SIO), and it is therefore non 
chaotic, confirming the role of inter-layer coupling in inducing (see also the Discussion section). 



Discussion 

The structure of the model 
Multi-layer architecture 

The reduction of the full multi-layer structure of primary visual cortex (a cartoon of which 
is shown in Figure [i]A) to a simpler two-layer network (Figure [Tj3) is a drastic simplification. 
Throughout this paper, we have emphasized that the two main cortical thalamo-recipient layers, 
i.e. IV and VI [42^|43j|45j are included within distinct model layers, corresponding respectively 
to the upper and the lower ring in our network architecture. We do not include separate rings 
for each of the six cortical layers. However, in order to reflect the poly-synaptic nature of the 

we 



pathway from cortical layer IV to VI — passing through layers II/III and V [45 46,49-51 
have made the latency of the connections from the upper to the lower model layer larger than 
for the connections from the lower to the upper model layer. The incorporation of additional 
layers within our model is in principle possible, but at the price of increasing further an already 
large number of parameters. Our choice of introducing just two layers was guided by the need 
to keep the model as simple as possible, while retaining a multi-layer structure. 

In the simulations described above, the external drive is smaller to the lower layer than 
to the upper layer. This choice was motivated by the fact that thalamic projections toward 
layer IV are more numerous than toward layer VI [45^. Nevertheless, it should be noted that 
layer VI neurons have dendritic arborizations extending into layer IV where they can receive 
additional thalamo-cortical inputs |69|. However, as illustrated in Appendix S8, the behavior of 
the network remains qualitatively the same, if one adopts identical external drives for the two 
layers. A second aspect that we have neglected about differences in the external drive to different 
layers, is the fact that the size of receptive fields depends on laminar location. In particular 



the receptive fields of layer VI neurons can be larger than the ones of layer IV neurons [70 , 71 
However, a proper description of the stimulus-size dependence of the inputs would require as 
well to take into account horizontal interactions between different layer IV receptive fields fitting 
into a same larger layer VI receptive field, a modeling aspect that we hope to address in future 
investigations. 

Connectivity 

In our model intra-layer excitation is modulated more strongly with angular distance than 
intra-layer inhibition. However, the probability of inhibitory connections is larger than the 
probability of excitatory connection at any angular distance (Figure [TJii and Table [5] for details). 
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In addition, we choose conductance parameters such that individual inhibitory PSPs are stronger 



than excitatory PSPs 72 . Thus, intra-layer inhibition dominates intra-layer excitation at any 
distance. As a consequence, in the regimes explored in this paper, recurrent interactions are not 
sufficient to generate a tuned response by themselves. However they sharpen the tuning already 
present in the spatially-patterned feed-forward LGN input. We use in the model probabilities of 



connection compatible with the wide ranges reported by 73,74 . Other studies, like |72l, find a 



larger probability of inhibitory connection. We verified however that the qualitative properties 
of the induced regimes of activity are preserved when inhibitory connections are consistently 
densified (see Appendix S8). 

The dominantly inhibitory nature of mutual local interactions is essential in our model for 
the emergence of prominent collective oscillatory behaviors in our network. Oscillations are 
generated by mutual delayed interactions between inhibitory neurons, according to a standard 



mechanism already described in 25,27-30 . In our model, excitatory neurons are not required 
for the generation of oscillations. Excitatory neurons are entrained by the oscillation paced by 
inhibitory cells. Indeed, if the activity of excitatory neurons is completely suppressed, or if 
synapses from excitatory to inhibitory neurons are removed, while increasing the drive to in- 
hibitory neurons in order to maintain their rate of activity unchanged, the oscillations continue to 
exist and their frequency increases of less than five percent (see Appendix S8). We mention here 
that an alternative scenario exists in which the inhibitory-to-excitatory-to-inhibitory neurons 



feedback loop plays an active role in the generation of synchronous oscillations 27,30,75-77 



In this scenario delayed inhibitory feedback is still the cause of the oscillation, but the delay 
arise from the disynaptic nature of effective mutually inhibitory interactions, leading to a slower 
collective frequency. However, the analysis conducted in Appendix S8 clarifies that the scenario 
implemented in our model relies primarily on inhibitory interneurons alone. 

Inter-layer connections in our model are as dense as intra-layer connections, but inter- layer 
excitation is more sharply modulated than intra-layer excitation. This results in a smooth 
arrangement of vertical excitatory synapses reminiscent of the organization of cortex into a 
continuum of anatomical columns without rigid boundaries |78|. This arrangement is critical 
for the fast temporal decorrelation of induced oscillations at high contrast (see below). 

Whereas the net inter-layer coupling is moderately excitatory in a local center, it is inhibitory 
in the surround, as a combined effect of the broad profile of inter-layer inhibition and of the 
fact that lower-to-upper excitation toward inhibitory neurons (i.e. disynaptic inhibition) is 
less sharply modulated than lower-to-upper layer excitation toward excitatory neurons. This is 
required in our model to account for the increase in mean firing rate observed in layer inactivation 
experiments |79| (case F = in our model). 



The low and high contrast regimes 

Most of the simulations described above were performed in networks with a significantly smaller 
number of neurons {Ne = 4000 excitatory neurons and Ni = 1000 inhibitory neurons per layer) 
than in a real hypercolumn in VI. However, we checked that our results are robust against 
increases in network size. In particular, this is the case for the existence of two dynamical 
regimes induced respectively by low and high contrast stimulations and for the two distinct 
mechanisms underlying the fast temporal decorrelation and broad-band spectral modulations in 
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these two regimes. 

In the low contrast regime, the dynamics are asynchronous. However, the network tends 
to resonate at a specific frequency, producing an increase of power in the gamma frequency 
band, without developing stable oscillations. Weakly coherent oscillatory modes are excited 
only transiently by local noise and then quickly damped. 

On the other hand, in the high contrast regime the network activity is synchronous. However 
the collective rhythm undergoes random variations in the time interval between consecutive 
activity episodes in the network. This temporal irregularity is not due to local noise (note 
that, in our model, recurrent inputs dominate over feed-forward inputs at low as well as at 
full contrast). It is produced intrinsically by the dynamics by virtue of the interaction between 
distinct oscillating populations localized in the two subnetworks representing different depths in 
the cortical section. This results in rapid temporal decorrelation of the induced activity. 

The contrast at which the transition between these two regimes takes place depends on the 
strength of fluctuations in the background noise. For our choice of parameters, the transition 
occurs for C ~ 10%. However, as discussed in detail in Appendix S9, if the variance in the LGN 
input current is increased consistently without changing its mean value, the transition can occur 
for an external drive, which is so large that it cannot be reached even for stimuli at full contrast. 
In such a condition, the induced activity is still asynchronous at high contrast and only transient 
oscillations can be detected, as in the recent modeling study by Mazzoni et al. f37^. 

It has been observed experimentally that the gamma-band synchronization of membrane 
potential fluctuations of nearby cells in VI is larger in visually-induced activity than in spon- 
taneous activity. Furthermore it is sustained over long stimulation durations, independently 
from stimulus properties or from the simultaneous observation of synchronized spiking activity. 
This leads to voltage crosscorrelograms with a manifest oscillatory component at gamma-range 
frequencies, damped quickly within only two or three oscillation cycles |80j. These observations 
are compatible with the occurrence of a transition between an asynchronous low contrast regime 
and a synchronous high contrast regime. Indeed, pairwise CCos between membrane potentials 
are small in the low contrast regime (Figures [3p and [5^), but large in the high contrast regime 
(Figure [4t) andlst)), even if spike CCos are always small, in agreement with many experimental 
reports [5, 13 81f[86| . We remark that if the dynamics at high contrast would be asynchronous 
as the dynamics in absence of stimuli or for low contrast stimuli, then the pairwise crosscorrela- 
tions of both spikes and voltages should be weak. Therefore, the coexistence of weak correlations 
between spikes with stronger correlations between membrane potentials (displaying furthermore 
a damping oscillatory component) is suggestive of the existence of a synchronous, rather than 
of an asynchronous, regime. The dynamics at high contrast of our model, characterized by 
irregular spiking (leading to weak spike crosscorrelations) and by temporally irregular collec- 
tive oscillations (leading to quickly damped oscillatory voltage crossocorrelograms) is therefore 



compatible qualitatively with the experimental regime observed in 80 . Conversely, this com- 
patibility could not be claimed for the other two types of induced dynamics that our model can 
generate at full contrast, i.e. asynchronous, in the case of a large variance noise, or synchronous 
but approximately periodic (and therefore too slowly decorrelating) , in the case of suppressed 
inter-layer interactions (F = 0). 
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Synchronous chaos underhes the temporal decorrelation of the network col- 
lective oscillations in the high contrast regime 

The rapid loss of temporal coherence of the synchronous induced activity at high contrast is a 
remarkable property of our model. Features of the model such as inter-layer inhibition, asym- 
metric interaction latencies in the lower-to-upper or in the upper-to-lower direction or different 
LGN driving levels to the different layers are not required for this decoherence to occur. In con- 
trast, the strong local inhibition responsible for the local generation of the rhythm within each 
layer and the net excitatory interactions between neurons in close vertical alignment are crucial 
for this to occur. In fact, if the inter-layer excitation profile is altered by suppressing its mod- 
ulation with orientation distance while keeping its average strength constant, the decorrelation 
does not take place (see Appendix S8). 

A similar mechanism underlies the temporal decorrelation of synchronous oscillations in the 



network models studied by 39-41 . These papers showed that collective oscillations induced in 
two populations of neurons by local delayed inhibitory feedback can lose coherence when the 
two populations interact in an excitatory manner. In [4lj, we studied a rate model consisting 
of two networks, each composed of one excitatory and one inhibitory populations. Each of the 
networks was able to sustain synchronous oscillatory activity by virtue of the local inhibition. 
We computed the maximum Lyapunov exponent of the system (see e.g. [66]) to show that 
it undergoes a transition to a chaotic dynamical state when the two networks are coupled 
by sufficiently strong excitatory connections. In this state the network displays synchronous 
activity, but instead of being periodic, the temporal variations of the network activity are chaotic 
and thus the oscillations that the network tends to develop lose temporal coherence within a 
few cycles. A network operating in such a regime is said to be in a synchronous chaotic state. 



In [39 40 a single ring network with strong local inhibition was considered. The decoherence 
of the oscillations occurred as the network underwent a spontaneous clustering into groups of 
oscillating neurons effectively interacting in an excitatory manner. 

In agreement with the positivity of its largest Lyapunov exponent, also the dynamics of our 
hypercolumn model in the high contrast regime displays typical features of chaos: exponen- 
tially fast damping of the local oscillations autocorrelograms (Figure j7p,D), spreading of the 
oscillation-related power over an extended continuous interval (Figure and extreme sensitiv- 



ity to initial conditions (Figure [lOp and Figure 11 ). Therefore the decoherence of the population 



activity which occurs at high contrast stems in the present model from the fact that the network 
operates in a synchronous chaotic regime. We cannot exclude, obviously, that other mechanisms 
are contributing to the decorrelation of synchronous cortical oscillations. We stress nevertheless 
that such a global decorrelation, characterized by the coexistence of elevated instantaneous syn- 
chrony and fast loss of collective phase autocoherence, could not be induced by local external 
noisy inputs, unless they are spatially correlated over a range matching the size of the local 
circuit which generates the ongoing soscillation. 

We also conjecture that the underlying mechanism of synchronous chaos is very general 
as it occurs in models in which neurons are described in term of rate, integrate-and-fire or 
conductance-based dynamics, with a simplified as well as more complex multi-layer network 
architecture. We also conjecture that a similar mechanism should act in even more realistic 
models, incorporating for instance a two-dimensional spatial structure, similarly to the one used 
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87 , 88 , provided that local inhibition is strong enough to induce local oscillations and that 



excitation couples these local oscillators at a longer range. 
Comparison with previous works 

Chaotic dynamics as well as stable chaotic-like dynamics can occur in asynchronous states of 
activity |89| - |97| . In this cases, the network dynamics explores a high-dimensional manifold in 
the phase-space, while, in our model, the irregular sparse firing of many neurons give rise to 



collective synchronous chaos (SC) with a lower dimensionality 98 - 100| (the fractal dimension 
of the chaotic attractor is likely to be smaller than five, as discussed in the Appendix SIO). 

SC has also been found in previous models of local circuits in VI which consisted of only 
one single network with a ring architecture. The model studied by Hansel and Sompolinsky 
in [101 considered one neuronal population coupled with excitatory instantaneous synapses. 



It displayed a SC state in some appropriate range of parameters. However, in this model, SC 
was sensitive to the incorporation of synaptic time constants since it was destroyed with the 
introduction of synaptic time constants as small as 0.5 ms. The model by the same authors 
considered in [34^ considered two populations of neurons, one excitatory and one inhibitory, 
coupled via synapses with realistic synaptic time constants. The dynamics of the neurons were 
based on a Hodgkin-Huxley type model with several cellular and synaptic conductances. The 
pattern of connectivity had a "Mexican hat" with local excitation and broad range inhibition. 
Numerical simulations of the model showed that in an appropriate parameter range, the network 
settled in a SC state, characterized by strong temporal variability of the neural activity which 
was correlated across the hypercolumn. 

In both of these models, the SC state was characterized by strong neuronal pairwise spike 
correlations and wide variability in the firing of individual neurons which was induced by the 
chaotic nature of the population activity. This is essentially different from what happens in 
our two layers hypercolumn, in which, in the SC state at high contrast, the spike pairwise 
correlations are only slightly larger than in the low contrast asynchronous state, whereas the 
degree of irregularity in the spike trains are similarly large in both states {CV ~ 0.9). As a 
matter of fact, in the present model, the spike train irregularities are mostly due to the local 
noise generated by the external inputs and to a lesser extent by the internal dynamics. Voltage 
CCos are large due to macroscopically correlated chaotic sub-threshold fluctuations, but spike 



CCos are still small. Another essential difference is that in 34 the excitation was local and 
inhibition was broad, whereas the opposite is required in the present model, as well as in the 
single ring model in ^40j. Last but not least, it is not clear to what extent the chaotic dynamics 



found in 34, lOll were specific to the model adopted there for the single neuron dynamics. 



Predictions and perspectives 

The increase in synchrony of the activity with the contrast displayed by our model is in agreement 
with experimental results reported recently in monkey VI [9, 18 . More generally we should 



expect that varying a feature of a stimulus in a way that increases the external drive on VI 
network should have a similar effect. This is consistent with other recent results showing that 
varying the size of a visual stimulus |12| or attention [6, 11 strengthens the coherence in the 
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activity of VI neurons. 

In the low and large contrast regimes identified in our model the increased gamma power 
in the LFP spectra is broadband. At low contrast, the loss of coherence of the oscillations in 
the LFP in a few tenths of a milliseconds is due to noise. At large contrast, it is a consequence 
of the chaoticity of the LFP time-series. The behavior of our model in both these regimes is 



compatible with recent results by Burns et al 63 , because of its lack of sustained auto-coherence 
of induced oscillations. 

Our simulations predict that infra-granular layer inactivation should globally affect the exper- 
imentally observed spectral properties of induced LFPs by enhancing its periodicity. Single-layer 

or 



inactivation experiments based on pharmacological or local cooling techniques 79,102,103 
with optogenetic techniques |'23','24l might be used to test this prediction. Furthermore, manip- 
ulations in which the firing of a single additional spike is induced (or suppressed, analogously to 



the simulation of Figure 11 ) can be performed. Extreme sensitivity to single-spike perturbations 



was experimentally proved in the case of asynchronous spontaneous cortical dynamics [104]. It 
would be interesting to repeat similar experiments in a stimulus-induced regime of oscillatory 
activity, in order to study the impact of the addition of a single spike on the time-course of 
ongoing LFPs. 

In the present study we focused on the role of the interactions between cortical layers in 
promoting temporal decoherence of gamma oscillations via the generation of synchronous chaos 
in a network with the size of a typical classical receptive field in VI. It would be interesting 
to investigate whether horizontal interactions which extend at distances beyond the classical 
receptive field also contribute to the loss of temporal coherence via a similar mechanism when 
the visual stimuli are extended. The basic two-ring network developed in this paper can be 
replicated into a bi-dimensional architecture including long-range excitatory interactions in order 
to investigate this potential contribution. This framework can be also applied to assess how the 
phase relationship between activity at different locations in VI (e.g. between center and surround 
of an extended stimulus) depend on the polarity of long range interactions. Furthermore, an 
additional source of decorrelation might be inter-areal interactions occurring at an even longer 
range. 

Finally, we have here considered temporal decorrelation induced by excitatory interactions 
between populations oscillating due to delayed mutual inhibition. It would be interesting to 
investigate whether a similar decorrelation phenomenon can arise when the mechanism for the 
local generation of oscillations is different, and is based for instance on circuit loops with active 



involvement of pyramidal cells 27 , 30 , 75 - 77 , 105 , 106 



Methods 

Our model of a functional hypercolumn in VI consists of two interacting rings of neurons, an 
upper and a lower ring, each comprising Ne excitatory and Ni inhibitory neurons connected 
recurrently. We denote by A'^ = 2(Ne + Nj) the total number of neurons in the network. Each 
neuron is labeled by its location on the ring to which it belongs; i.e. by an angular coordinate 



•d, ranging conventionally from -90 to +90 degrees 53 ,541. All the neurons receive an external 



input composed of two contributions. One represents the LGN input to VI. It depends on two 
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parameters C and Ostim corresponding to the contrast and the orientation of a visual stimulus. 
The other contribution accounts for the background inputs VI receives from subcortical regions. 



Single neuron dynamics 

Throughout the paper, we use single-compartment Exponential Integrate-and-Fire model neu- 
rons (EIF; [lOTj). In this model the membrane potential V is given by the equation: 



at Tm C 



where C is the membrane capacitance, the membrane time-constant, Vl the leak potential, 
Isyn the total synaptic input current to the neuron. The function tpi^) is: 

At fV-Vr 

ipiV) = exp — 

For a constant input above a threshold current (~ 0.113 nA for the parameters adopted here) 
the solution of ([T]) diverges to infinity in finite time. This divergence is identified with the firing 
of a spike. The parameters and Vt characterize how sharp the initiation of the spike is and 
the voltage at which it occurs. The spike downswing is not explicitly modeled. After each spike 
event, the voltage needs to be reset. A refractory period must then follow. 

We model this refractoriness in a different way for excitatory and inhibitory neurons. In the 
case of excitatory neurons, following the emission of a spike at time t spikes the parameters Tm, 
Vt and Vl are updated according to the equations 1108 , 



Vl = vl + Av^ exp 



t tspike 



Bvl exp 



t tspike 



Vt = Vt + Avj. exp 



t - 1 



spike 



TA,Vt 



1 



1 



+ exp 



t - 1 



spike 



(3) 



(4) 



(5) 



The membrane potential is reset to a value Vreset which is sub-threshold. Furthermore Vt is 
strongly depolarized after a spike. Therefore the event that two spikes are closely emitted in 
time by a same neuron is extremely unlikely and, in practice, never occurs. 

For inhibitory interneurons, we use a "hard" refractory period instead, suspending the nu- 
merical integration for a time r^e/ after voltage reset 
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Therefore, Vl 



Vl, Vt 



Vl and 



Parameters for excitatory neurons are chosen to coincide with fits of pyramidal neurons 
traces, following |108 . We use analogous parameters for inhibitory neurons, apart from halved 
membrane capacitance and time constant Tm, consistent with experimental evidence 109] and 



fits of interneuronal traces presented in 110| . All single neuron parameters are given in Tables [l] 
and [3 
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The synapses 



We use three kinds of synaptic currents, modeling inhibitory (GABA-type), fast excitatory 
(AMPA-type) and slow excitatory (NMDA-type) synaptic inputs. No voltage dependence is 
introduced for the parameters of the slow excitatory synaptic current. A spike in an inhibitory 
pre-synaptic neuron evokes a GABA-type post-synaptic potential (PSP) in all the post-synaptic 
neurons; a spike in an excitatory presynaptic neurons evokes composite AMPA- and NMDA-type 
PSPs. 

The synaptic current produced by a single incoming spike is described as Isyn,spikei't) = 
—dsyniV — Vsyn)s{t), where gsyn is the peak synaptic conductance, Vgyn the reversal potential 
of the synapse {Vampa = Vnmda = 0.0 mV, Vqaba = —75 mV). Denoting as ispjfce the time 
of pre-synaptic firing and with d the synaptic latency, the function s{t) is: 



s{t) 



1 



'~(*spifce + rf) 

e — e~ 



t-(t. 



pike 



+ d) 



(6) 



where the constant N is such that it normalizes to unity the peak of s{t). All the synaptic 
conductances in the network are calibrated to give unitary PSPs at resting potential in a range 



compatible with experimental observations 72 



The values of the synaptic times and synaptic peak conductances are given in Table [3], for a 
network including A'^^; = 4000 excitatory neurons and Nj = 1000 inhibitory neurons per layer. 
Synaptic peak conductances are rescaled for larger networks, according to Eqs. ([sj) and ([9|. 
Synaptic latencies are given in Table |4} 



Network connectivity 



Each of the two layers of the hypercolumn is modeled by a ring -network (34|[53|[54|[lll]. Un- 
less specified otherwise, the simulations described in this paper were performed for a network 
comprising Ne = 4000 excitatory cells and Ni = 1000 inhibitory cells per ring, for a total 
of = 2{Ne + Ni) = 10000 neurons in the hypercolumn. Note that a very similar network 



architecture was used in [112 , 1131 but with a completely different interpretation. 



Intra-layer and inter-layer excitatory and inhibitory connections are random. The proba- 
bility of connection between two neurons is spatially modulated and depends on the angular 
coordinates 'dpre and "dpost of the pre- and post-synaptic neurons. It also depends on the nature 
(excitatory or inhibitory) of pre- and post-synaptic cells and on their absolute (lower or upper 
layer) and relative (intra-layer or inter-layer) depth. All the profiles of connection probability 
are parameterized as: 



) "post ) 



pre 



'&post ) 



(7) 



Here, [•]+ denotes rectification; i.e. = if x < 0, else x]+ = x. The probabilities of 

connection for intra-layer excitatory and inhibitory connection are identical for each of the two 
layers. 

In order to study the scaling properties of the dynamics it is important to guarantee that the 
spatial mean and spatial fluctuations of the time averaged recurrent synaptic currents received by 
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each neuron are preserved when considering networks of different sizes. This requires a suitable 
modification of the probabihties of connection and of the peak synaptic conductances when 
passing from a network of size to a network of size N' [35]. For an arbitrary peak recurrent 
synaptic conductance gx, the probabilities of connection (and, correlatively the average number 
of pre-synaptic cells of each type) are scaled as: 



and peak conductances as: 



PxNg, = P'.N'g', (9) 



Here the index x stands for different kinds of synaptic connections, each one potentially charac- 
terized by different mean probabilities of connections and connection strengths (i.e. originating 
from upper or lower layer excitatory or inhibitory neurons and directed toward upper or lower 
layer excitatory and inhibitory neurons). 

Sizes between Ne = 1000 — for a total network size of = 3000 neurons — and A''^; = 32000 
— for a total network size of = 80000 neurons — are compared in scaling analysis of synchrony 
properties. The parameters for Ne = 4000 and Nj = 1000 are given in Table [5| Probabilities 



of connection are compatible with the ranges reported by 73 , 74 . 



Model of the LGN input 

We assume that the firing rate of a single LGN neuron, r(C) is related to the stimulus contrast, 
C, (C = 1 ^ 100 % ) by the equation [TTl] : 



r(C7) = ro + rilogio(l + C) Hz (10) 

where tq is the spontaneous activity of the neuron in dark conditions. Subsequently, we model 
the LGN input to a cortical cell as an AMPA-type synaptic connection with peak conductance 
9LGN, driven by homogeneous Poisson spike trains with rate RLGN{'&,'&stim,C), 

RLGN{^,^sUm,C) =Ro + [Rl{C) (1 - e + eCOs2(l? - l}sUm))]+ (11) 

with: 

i?i(C) = i?ilogio(l + C) (12) 

Here the parameter e controls the broadness of tuning of the LGN input. It is set to 1 in all our 
simulations. Note that Rlgn is maximum when 'dsUm = '&■ The contrast C and, correspondingly, 
the term Ri (C) can also be time-dependent (see later section on peristimulus time histograms) . 
The LGN input targets both layers. There is anatomical evidence that thalamo-cortical synapses 



target mainly layer IV and to a lesser extent layer VI 42 , 45 . Accordingly, in all the simulations 
presented in this paper, qlgn in the lower layer is smaller by a factor of two than in the upper 
layer. Parameters describing LGN input properties are given in Table [6} 

For the adopted parameters, feed-forward inputs from LGN never dominate over recurrent 
inputs from the two layers of the network, consistently with the massively larger number of 
cortico-cortical synapses than thalamo-cortical synapses in the primary visual cortex (45|. The 
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relative weight of feed- forward inputs with respect to recurrent inputs depends on contrast, 
doubling in our model from no more than 20% at low contrast to no more than 40% at high 
contrast stimulation (not shown). 

An alternative parameter choice for the tuned component of the LGN input, leading to 
noisy input current with a larger variance, is analyzed in Appendix S9. For the noisy inputs 
used in this paper as well as for the ones used in Appendix S9, the resulting sub-threshold voltage 
fluctuations are on the order of ~6 mV at full contrast, compatible with experimentally observed 
fluctuation ranges 114 115| . Voltage fluctuations are comparable in the two regimes, because 
the increase in amplitude of external input current fluctuations is paralleled by a decrease in 
amplitude of net input conductance fluctuations, due to reduced synchrony among the recurrent 
inputs (see Appendix S9). 

More details about the mapping from stimulus contrast to input rates can be found in 
Appendix Sll. 



Background cortical noise 

In addition to the LGN input, excitatory and inhibitory cells are driven by an untuned noisy 
input, representing the background firing of other cortical areas. This input is modeled by a 
single AMPA-type synapse per cell, with peak conductance gbg activated by Ornstein-Uhlenbeck 
processes [61] . Input spikes are generated independently for each cell; however all the cells share 
the same instantaneous input rate Rbg{t) obeying the stochastic differential equation: 

= --(R^At) - /^..) + \f^m (13) 

at Tbg y Thg 

where ^{t) stands for Gaussian white noise and fx^g is the mean, ai,g the volatility and r^g the 
filtering time-constant of the stochastic process. Parameters are given in Table [Tj 



Numerical integration scheme 

The dynamical equations are integrated using a fourth-order non-adaptive step Runge-Kutta 
scheme. Integration step was 0.2 ms. Because of the exponentially fast divergence of the 
membrane in relation with firing, particular care is needed to ensure the stability of the numerical 
integration of equation ([T]) . Following [107j , the numerical integration of the membrane potential 
y of a given neuron is stopped as soon as V reaches a finite cutoff voltage V*. In our simulations, 
we use V* = —30 mV. This choice ensures that the non-linear term ip{V*) is the dominant 
contribution to the neuronal currents for V > Vth- Under this condition, the leakage and the 
synaptic currents can be neglected, making it possible to compute analytically the time left before 
the actual divergence of the potential. Assuming that the integration is stopped at t = tgtop 
when V = Vstop > V*, the time of the next spike is given by t^pifce — tstop + Tm,e'^^^~^''*°p)/^^. 
In addition, for our choice of V*, t spike — tstop is large compared to the integration-step At, thus 
avoiding numerical errors in spike-time estimation due to the exponentially fast growth of V in 
proximity of the divergence. The membrane potential is then reset to a value Vreset immediately 
after a spike. 
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The Ornstein-Uhlenbeck process giving Rbg{t) is computed using the properties of the exact 
solution to equation ( (131 ). This means that Rbg{t + At) is normahy distributed wi th m ean 
Rhg{t)e-^^/^b9 + ^^^(1 _ e-^*/^''9) and standard deviation s = abg^/l- e-^^^/^f-g 



m 
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Response tuning and contrast response function 

In order to study the tuning properties of the neuronal responses we present stimuli at 12 
different orientations 'dsUm in an interval between -90 degrees and +90 degrees, at at least five 
different contrast values C per each orientation. 

Tuning curves are derived for each neuron by measuring their average firing rate for each 
of the tested orientations and contrasts and are characterized by computing their skewness and 
their circular variance |117| . See Appendix S3 for more details. Population average tuning 
curves are computed after rotating single neuron tuning curves so that their maximum is at 
■(? = (see Figure [2]). 

The contrast-response functions, CRF(C), are computed for each neuron by measuring its 
peak firing rate (i.e. its firing response to a preferred orientation stimulus) at each given level 
of contrasts. Each individual CRF is fitted to a hyperbolic ratio function (60|: 



CRF(C) = i?^,,— -— (14) 

Measures of synchrony 

To measure the degree of macroscopic synchrony in the steady state of a network comprising 



an arbitrary number of neurons, we follow the method used in 34,35 . It is grounded on 
analysis of the temporal fluctuations of macroscopic observables of the networks such as the 
instantaneous activity or the instantaneous membrane potential averaged over a population of 
neurons of size K. For instance, for the latter, one evaluates at a given time, t, the quantity: 

K 

K 



The variance of the time fluctuations of V{t) is 

^v = {[V{t)f)-[{V{t))f (16) 



where (...) = dt . . . denotes time-averaging over a large time, T. After normalization of ay 
to the average over the population of the single cell membrane potentials: 

^v. = (mt)?) - mm' (i?) 

one defines a synchrony measure, x{K) by: 

X\K) = ^^^ (18) 
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This measure takes values between and 1. In the Hmit — )• oo it behaves as: 

xW = x(oo) + ^ + 0(l) (19) 



where a is some constant, between and 1. In particular, x(^) = li if the system is fully 
synchronized (i.e., Vi(t) = Vii) for all i), and X (^) = 0(1/y^K)) if the state of the system 
is asynchronous. Asynchronous and synchronous states are unambiguously characterized in the 
thermodynamic limit (i.e., when the number of neurons is infinite). In the asynchronous state, 
x(oo) = 0. By contrast, in synchronous states, x(oo) > 0- 

To characterize the degree of synchrony in the membrane potentials of neurons i and j, we 
compute the cross-correlation function: 

The value of the normalized cross-correlogram for zero time-lag gives the pairwise crosscorrela- 
tion coefficient (CCo): 

cco(yi,y,) = cc(yi,yj)[T = o] (21) 

To estimate the degree of synchrony in the spiking activity of these neurons, discrete spike trains 
are first convolved with a square window of width i?, thus generating a continuous spike-count 



signal. Equations (20) and (20), with Vi replaced by such smoothed spike trains, is used to 



compute crosscorrelograms and CCos for spiking activities 118| . We use a smoothing window 
size of S = 20 ms. 

CCos and crosscorrelograms are estimated over simulated recordings lasting 100 s of real time. 
For CCos between membrane potentials only pairs of neurons within a 18° region centered on 
an angular coordinate matching the orientation of the presented stimulus are considered. In the 
case of spike trains, neurons in this region whose spike train contained fewer than 100 spikes are 
further excluded. Various stimulus orientations are pooled together to improve the estimation 
of CCo distributions. 

Local field potentials 

LFPs are believed to be an aggregate measure of the synaptic activity of several hundreds of 
neurons in the vicinity of the recording electrode [119 , 120 . To evaluate the LFP in a given site, 



we thus arbitrarily average input synaptic currents in a small angular sector of 9° centered on the 
considered angular position. LFPs are estimated over neurons of the upper layer only, reflecting 
the fact that superficial neurons should supply the largest contribution to the signal recorded 
by an applied recording tip. For the normally used size of iV^ = 4000 excitatory neurons and 
Ni = 1000 inhibitory neurons per layer, this corresponds to averaging over 200 excitatory and 
50 inhibitory upper layer neurons for each considered LFP recording site. 
Autocorrelograms of the LFPs are computed as: 



AC[LFP](t) = (LFP(t) • LFP(t + r)) - (LFP(i))2 



(22) 
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We evaluate non-normalized autocorrelograms, in such a way that the zero-lag value AC[LFP](r = 
0) measures the variance of the temporal fluctuations of the LFP and has known size-scaling 
properties, which are different in synchronous and asynchronous regimes [35j. 

Power spectra are computed using conventional FFT techniques, as the square modulus of 
the Fourier Transform of signal autocorrelation. Windowing is applied to LFP-like signal time- 



series to reduce unwanted frequency leakage, following the Welch method 121| ). An additional 
moving average smoothing is applied for visualization purposes. We measure power in arbitrary 
logarithmic units. Since we are interested in qualitative analysis of the overall shape of the 
spectra rather than in absolute power estimations, for each considered regime we assign a unit 
value at the power at Hz frequency for 0% of contrast. 

Autocorrelation and spectral analysis of LFP-like signals are based on time-series lasting 100 
s of real time, with a sampling rate of 5 kHz. 

Multi-unit-activity 

The MUA signal reflects the spiking activity of few neurons in the immediate vicinity of the 



recording electrode 122 . Typically, the recorded MUA is separated in only a small number of 



contributing single units 123 . To evaluate MUA at a given site, following fTT3|, we sum together 
the spike trains of three randomly selected cells within a small angular sector of 9° centered on 
the considered angular position (the same used for the evalutation of the LFP). This discrete 
signal is then transformed into a continuous signal by convolving it with a gaussian window (1 
ms of variance). 



We compute then the coherence 124 between the LFP and the MUA at a same site by 



taking the modulus of the normalized product of their complex Fourier Transform, using the 



Welch method 1211, as in the case of the LFP power spectrum estimation: 



SlfpU) if) 

V\SLFp{fr-\SMUA(.W 

where SlfpH) and SMUA^f) denote the Fourier transform of the autocorrelograms of the LFP 
and of the MUA signals, respectively, a star denotes complex conjugation and | • | complex 
absolute value. Such MUA-LFP coherence is a real quantity in the unit interval < C(/) < 1, 
and provides an absolute (linear) measure of the phase synchronization between the two signals 
in different frequency bands. We average then the result over twenty different randomly chosen 
triplets of cells, in agreement with the experimental habit to average together different MUA 



recordings with only approximately similar selectivity properties 1113 



Note that both our MUA and our LFP signals are not designed to reproduce faithfully 
realistic electrophysiological signals. Due to the extreme sensitivity of the coherence measure 
[124 , therefore, we cannot expect the MUA-LFP coherence measure we compute to have nothing 



more than an illustrative purpose. 

Inter-layer coupling strength and layer decoupling 

Layer decoupling experiments are performed by multiplying the peak conductances of all the 
AMPA-type, NMDA-type and GABA-type synapses from the lower layer to the upper layer by 
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a factor T varying between 1 and 0. A value of 1 corresponds to the case of fully- interacting 
layers, and a value of corresponds to fully uncoupled layers. 

For each excitatory neuron in the upper layer, at high contrast, the peak response after 
layer decoupling is compared with the peak response of the same neuron in the fully coupled 
network case. In comparing peak responses, we incorporate the fact that the tuning curves of 
many neurons change their preferred orientation or their skewness after full or partial lower layer 
inactivation. 



Peristimulus time histograms 

To simulate the flashing of a grating, for a given network realization we perform numerical 
simulations in which the tuned LGN input rate is not constant. More specifically, this tuned 



component is still modeled according to equations (11) and (12), but the contrast is now mod- 
ulated in time: 

^_^m_/0% n-1.5s<t< (n-1.5 + 0.5) s r> - \ 2 (24) 

^ -^W- I 95% (^.1.5 + 0.5) s<t< (n + l)-1.5s (^4j 

Phases lasting 0.5 s in which Ri = 0.0 are therefore alternated with phases lasting 1 s in which 
Ri = i?i(95%), leading to a square wave time-course of the input LGN rate. We consider only 
cells whose preferred orientation falls within a sector 9° wide centered on the orientation of the 
presented stimulus and we use four different stimulus orientations. For each of the orientations, 
the stimulus is flashed 1000 times. An overall sample of 800 cells (200 per orientation) is thus 
considered. We count across stimulus repetitions and cells how many times a spike is emitted 
within 2 ms from a specified time following the onset of a stimulus. The probability that a neuron 
will fire at a given time after stimulus presentation is then evaluated as (number of spikes in a 
time-from-stimulus bin) / (number of stimuli repetitions) / (number of sampled cells). 



Single spike perturbation 

To study the sensitivity of induced dynamics to a small perturbation, we perform a simulation in 
which just a single spiking even is omitted and we compare it with the unperturbed simulation. 
We select a putative spiking time of a neuron whose preferred orientation matches the one of 
the applied 95% of contrast stimulus. No spike is then sent to the post-synaptic targets, we only 
reset the potential and the other time-dependent parameters of the failing presynaptic neuron 
to their just-after-spike values. Precisely the same realizations for all the stochastic noisy input 
processes are taken for the unperturbed and the perturbed dynamics. 



Estimation of the largest Lyapunov exponent 

Rather than estimating the largest Lyapunov exponent through a "direct" method based on the 
application to the system of vanishingly small perturbations or through a "exact" method based 



on the integration of the linearized system 125 , we measure the largest Lyapunov exponent 



of the induced dynamics of the system at high contrast through a non-linear analysis of a long 
time-series (600 minutes of real time) of the associated LFP signal. For a thorough introduction 
and a rationale to the used methodologies the reader is invited to refer to textbooks like |68j . 
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We try just here to give a flavor of the employed techniques. The first step is the construction 
of proper "embeddings" of this time-series. Given a discretely sampled time-series it = LFP(t), 
a reconstruction delay r and an embedding dimension m, we construct a new m-dimensional 
sequence: 

f-t = (4-(m-l)r> h~(m-2)T: ■ ■ ■ At) (25) 



It can be proven 126, 1271 that the latter time-delay embedding provide in general a one-to-one 



image of the original phase-space attractor of the dynamics generating the measured time-series, 
provided that the used embedding dimension m is large enough. The general idea of the method 
is then to identify by systematic search pairs of points it and if which lie at a (euclidean) 
distance in the delay-embedding space smaller than a specified very small e. Such points are 
said to be neighbors. It is therefore possible to consider the distance (5o = — ^t' || as a "small 
perturbation" , which should grow exponentially in time if the dynamics is chaotic. The eventual 
divergence of the trajectories originating by neighbor points can be monitored by the quantity 
5k = \\it+k — it'+k\\- If a time range for which 5k oc (5oexp(AA;) than A cohincides with the 



maximum Lyapunov exponent Xmax 

128 129 



We select a reconstruction delay of r = 400 ms, much larger than the decorrelation time of 
the induced LFP oscillation. The minimum embedding dimension for a consistent estimation of 
the largest Lyapunov exponent can be estimated by monitoring the fraction of "false neighbors" 
pairs, i.e. pairs of points that are neighbors in a Z?- dimensional embedding (due to a projection 
of the attractor to a space with a too small dimensionality) but that there are no more such in 
an embedding with a larger dimension D' [130J . Such analysis, summarized in the Appendix 
SIO, indicates a critical embedding dimension lower than five and probably larger than three 
(even if a precise estimation is difficult due to the presence of noise). Practically, we estimate 
the largest Lyapunov exponent by evaluating the quantity: 

ln|(e,m,fc) = /ln(-L^ ll^m " | ) (26) 

\ V it' 6Wt(^) 

for various e, m and k, where Ut{e) is the set of points at a distance d < e from it and {■)t 
denotes average over time. If In |^(e, fn^ k) displays a linear increase in a reasonable range of k 
with matching slopes for different sufficiently large m and for few decades of e, than the average 
slope of the linear sections of In |^(e, tn, k) provides a robust estimation of Xmax- More details 
on our estimation of Xmax for the high contrast induced LFPs for P = 1 and P = are given in 
the Appendix SIO. 
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Figure captions 

Figure 1. Schematic drawing of the model hypercolumns. A: cartoon of the loop circuit 
among the 6 layers of striate cortex. Thalamo-recipient layers are indicated by pink shading. 
B: two-rings network, corresponding to a hypercolumn with interacting layers. LGN inputs 
are weaker toward the lower layer than toward the upper layer. C: the single ring network 
for each layer of the model hypercolumn. LGN inputs target both excitatory and inhibitory 
neurons. D: spatial profile of LGN input. E: spatial modulation of the probability of connections 
between two cells in the same layer, separated by an angular distance At?. Red line: excitatory 
connections. Blue line: inhibitory connections. F: spatial modulation of the probability of 
connections between two cells in different layers, separated by an angular distance A-d. Red line: 
upper-to-lower layer excitatory connections and lower-to-upper excitatory connections toward 
excitatory neurons. Magenta line: lower-to-upper layer excitatory connections toward inhibitory 
neurons. Blue line: lower-to-upper and upper-to-lower layer inhibitory connections. 

Figure 2. Response tuning and contrast response. A: tuning curves for different 
contrast levels (re-centered average over Ne = 4000 excitatory neurons in upper layer). Solid 
lines represent Gaussian fits. B: contrast response functions (blue curve: average over Nj = 1000 
inhibitory neurons in the upper layer; red curve: average over Ne = 4000 excitatory neurons in 
the upper layer). Solid lines represent hyperbolic ratio fits. 

Figure 3. Low contrast dynamics. Dynamics of the upper layer for the presentation 
of a 2%-contrast stimulus. A: raster plot of the excitatory population activity and associated 
time-histogram of the rate of spiking cells. The histogram bar heights denote the fraction of 
upper layer excitatory cells that fire in the bin. Bin-size is 2 ms. B: spike trains of 6 excitatory 
cells highly activated by the presented stimulus. C: membrane potential traces for two neurons 
stimulated simultaneously at close-to-preferred orientation (2 top neurons of Panel B in red 
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and green). D: pairwise correlations between spike trains (left, cyan histogram) and membrane 
potentials (right, blue histogram) of highly active neurons. 

Figure 4. High contrast dynamics. Dynamics of the upper layer for the presentation 
of a 95%-contrast stimulus. A: raster plot of the excitatory population activity and associated 
time- histogram of the rate of spiking cells. The histogram bar heights denote the fraction of 
upper layer excitatory cells that fire in the bin. Bin-size is 2 ms. B: spike trains of 6 excitatory 
cells highly activated by the presented stimulus. C: membrane potential traces for two neurons 
stimulated simultaneously at close-to-preferred orientation (2 top neurons of Panel B in red 
and green). D: pairwise correlations between spike trains (magenta histogram) and membrane 
potentials (red histogram) of highly active neurons. 

Figure 5. Pairwise crosscorrelations of spike trains and membrane potentials. 
Autocorrelograms and pairwise crosscorrelograms of spiking activity and membrane potentials 
for three upper layer excitatory neurons. A: spiking activity, low contrast, C = 2%. B: membrane 
potential, low contrast, C = 2%. C: spiking activity, high contrast, C = 95%. D: membrane 
potential, high contrast, C = 95%. Auto- and crosscorrelograms are normalized (for zero time- 
lag, autocorrelograms peak at one and crosscorrelograms at the correlation coefficient). The 
units for the time-lag axis are ms. Colors are as in Figures [3p and [4p. Rows and columns 
correspond to different neurons. The angular coordinates of the three neurons are 0°, —10° and 
10°. 

Figure 6. The Measure of synchrony as a function of the contrast and different 
network size. A: The synchrony measure, X) increases abruptly with the stimulus contrast 
N = 10000 (solid line) and N = 40000 dotted line). B: The synchrony measure x as a function 
of the network size for spontaneous activity (zero contrast, grey line), low contrast (blue line) 
and high contrast (red line). The dashed line corresponds to a power- law decay with exponent 
-0.5. 

Figure 7. The autocorrelograms of the local field potentials. A-B: low contrast, 
C = 2%. C-D: high contrast, C = 95%. Scalings of non-normalized autocorrelograms are 
shown in B and D. In both cases the damping of secondary peaks is faster for larger network 
sizes. Zero-lag autocorrelation vanishes for large sizes at low contrast but not at high contrasts. 
Non-normalized autocorrelations are measured in nA^. 

Figure 8. Spectral properties of the LFP and MUA for diff^erent contrasts. A: The 
power spectra are shown for the LFP induced by a stimulus at preferred orientation. Isolated 
peaks do not appear even for very high contrast stimuli. B: The average coherence spectra are 
shown between the MUA and the LFP induced at a same location by a stimulus at preferred 
orientation. MUA-LFP coherence and LFP power are modulated by contrast changes in the 
same broad frequency range in the gamma band (30-100 Hz). 

Figure 9. Effects of the layer decoupling on the dynamics of the hypercolumn. 
Changes for decreasing inter-layer coupling and for a stimulus at high contrast with preferred 
orientation. A: population average peak firing rate for the excitatory neurons in the upper layer. 
B: synchrony level x- C: autocorrelograms of LFPs for intermediate strengths of the inter-layer 
coupling (r = 0.8,0.6 and 0.2). D: corresponding LFP power spectra. E: autocorrelograms of 
LFP for preferred stimulation at high contrast for the case of fully uncoupled layers (F = 0.0). 
F: corresponding LFP power spectrum. Spectra are also plotted for lower levels of contrast and 
are characterized by a narrow peak at a contrast-dependent frequency. 



36 



Figure 10. Short-term response. Population firing responses to repeated presentations 
of a high contrast stimulus for fully coupled layers (A-B, F = 1) and for fully uncoupled layers 
(C-D, r = 0). A and C: peristimulus-time (PST) histograms, based on the firing responses of 
500 cells to 1000 presentations of stimuli with optimal (or close to optimal) orientation. B and 
D: examples of upper layer excitatory population responses for three presentations of the same 
stimulus. 

Figure 11. Chaotic sensitivity to a single spike perturbation. A black triangle de- 
notes the time of a small perturbation to the network dynamics (for 95% of contrast stimulus and 
fully-coupled layers, F = 1), in which a single spiking event is omitted. Already after the second 

oscillation cycle, the unperturbed and perturbed population dynamics have diverged, as visual- 
ized by the raster plot (A) and the population rate histogram (B) of the upper layer excitatory 
population. Blue color denotes unperturbed dynamics and red color perturbed dynamics. 

Supporting information: captions 

Supporting informations can be downloaded at the web address: 

http : //www . ploscompbiol . org/article/inf o7.3Adoi7.2F10 . 137iy.2Fjournal.pcbi . 1002176 

Figure SI. CV and firing rate distributions. Distributions of CVs and firing rates for highly 
active upper layer excitatory neurons (orientation preference within ±5° from stimulus) arc here 
shown for the spontaneous activity (C = 0%), for the low contrast regime (C = 2%) and for 
the high contrast regime (C = 95%). Distributions of CVs (A-C) and of firing rates (D-F), 
from top to bottom, for the spontaneous activity, for the low contrast regime and for the high 
contrast regime. 

Figure S2. Heterogeneity of single neuron tuning curves. A: re-centered single neuron 
tuning curves for 3 upper layer excitatory neurons. B: distribution of peak rates. C: distribution 
of tuning width. D: distribution of tuning skewness. Distributions are relative to the upper layer 
excitatory population. 

Figure S3. Contrast invariance of tuning width. Tuning curves normalized to peak height, 
for: fully coupled layers case (A), fully uncoupled layers case (B) and strong noise case (C; see 
Figures S13 and S14). In cases A and B, contrast invariance is only approximate and does not 
hold for weak contrasts. Contrast invariance at low contrasts is improved in the strong noise 
case, in agreement with theory. 

Figure S4. Heterogeneity of single neuron contrast response functions. A: single neuron 
CRFs for 3 upper layer excitatory neurons. B: distribution of saturation rates. C: distribution of 
CRF steepness. D: distribution of mid-range contrasts. Distributions are relative to the upper 
layer excitatory population. 

Figure S5. Correlation coefficients for spontaneous activity. A: pairwisc correlations (CCos) 
between membrane potentials. B: pairwise correlations between spike trains. CCos are computed 
as described in the Materials and Methods section of the main article. 
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Figure S6. Oscillatory structure of induced LFPs. A-B: example LFP traces from our 
model, evoked by a stimulus with 2% of contrast (A) or 95% of contrast (B). C-D: autocor- 
relograms and crosscorrelograms of evoked LFPs. For a fixed stimulus orientation, we monitor 
AGs and CCs of LFPs in regions responding preferentially to this stimulus orientation or to an 
orthogonal stimulus orientation. The analysis is performed for 2% of contrast stimuli (C) or 
for 95% of contrast stimuli (D). An oscillatory structure is present in LFP independently from 
spiking and is correlated over the entire ring. E-F: LFP temporal decorrelation at intermediate 
contrast levels (E: C = 20%; F: C = 4%). Damping of secondary peaks is fast at any contrast. 
Units are nA^. 

Figure S7. Dynamics of the lower layer. A: raster plots of the activity of the lower layer 
excitatory neurons (lower raster) and the upper layer excitatory neurons (upper raster) in the 
high contrast regime dynamics {C = 95%). B: the latency between induced oscillations in the 
upper and in the lower layer is estimated through the crosscorrelogram of LFPs in the two layers 
(high contrast regime, C = 95%). The upper layer advances the lower layer ~ 2.8 ms on average. 

Figure S8. Induced responses for fully uncoupled layers: changes in tuning curves. For P = 
(full layer uncoupling) , the dynamics of the upper layer is equivalent to the case where there is 
full inactivation of the lower layer. After layer uncoupling, and consistently with reference [79] 
(Allison and Bonds, 1994), we observe changes in preferred orientation, peak response rates and 
tuning curve width and skewness of single neuron tuning curves. We report here distributions 
of parameter changes (vertical dotted lines denote average parameters for fully coupled layers, 
P = 1). A: distribution of preferred orientation shifts. Preferred orientation of individual cells 
can move clockwise or anti-clockwise within a range of ~ ±30° but the distribution of shifts 
is symmetric, with no significant change at the population level. B: distribution of peak firing 
rate changes. The mean peak rate change is weakly positive, reflecting the overall inhibitory 
nature of inter-layer coupling. C: distribution of broadness changes. On average, the width of 
tuning curves is slightly increased. D: distribution of skewness changes. Skewness changes are 
observed in both directions and their distribution is symmetric, with no significant change at 
the population level. In general, the large heterogeneity in the effects of layer uncoupling on 
tuning properties must be noted. 

Figure S9. Induced responses for fully uncoupled layers: dynamical properties. Response of 
the upper layer for the presentation of a 95%-contrast stimulus. A: raster plot of the excitatory 
population activity and associated time-histogram of the rate of spiking cells. The histogram 
bar heights denote the fraction of upper layer excitatory cells which firing in the bin. Bin- 
size is 2 ms. B: spike trains of 6 excitatory cells highly activated by the presented stimulus. 
C: membrane potential traces for two neurons stimulated simultaneously at close-to-preferred 
orientation (2 top neurons of Panel B in red and green) . This dynamics is strongly synchronous 
and approximately periodic. For increasing network size, oscillations tend to become more 
periodic, and collective synchrony does not vanish (not shown). 

Figure SIO. Numerical experiments for chaos assessment. All the methods are described 
in Text SI. A: estimation of the minimal embedding dimension. The fraction of false neighbors 
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is plotted against the embedding dimension for a LFP time-series generated by a full contrast 
preferred orientation stimulus (F = 1, C = 95%). N = 1000 pairs of candidate neighbor points 
have been considered for each embedding dimension (e < 10~^). A threshold of R* = 10^ has 
been taken. A single LFP time-series long 10 hours of real time, with a sampling rate of 0.01 
ms has been used for the estimation. The resulting embedding dimension appears to be m > 4. 
B-C: extraction of the largest Lyapunov exponent \max for the dynamics induced by a full con- 
trast preferred orientation stimulus. The relative growth in time of the average separation 
between LFP trajectories originated from neighbor points is plotted against time, for various 
embedding dimensions (average over at least N = 1000 pairs of neighbors per considered em- 
bedding dimension). A section of exponentially fast growth (linear growth in a semilogarithmic 
plot, denoting deterministic chaos) is identified for sufficiently large embedding dimension in 
the case of a hypercolumn with interacting layers (F = 1, panel A), but not in the case of a 
hypercolumn without inter-layer interactions (F = 0, panel B). 

Figure SI 1. Alternative parameter choices: network with increased symmetry. We assumed 
in the main text that the LGN input to the lower layer is weaker. We show here results for 
the case in which the LGN input rate to lower layer is the same as to the upper layer and in 
which the latency of the upper-to-lower layer connections is as short as the latency of lower- 
to-upper layer connections. A: raster plot of the evoked activity of the upper layer excitatory 
population for a 95% level of contrast stimulus. B: autocorrelogram of the evoked LFP. Units 
are in nA^. Note that synchronous chaos is still present, as evidenced by the fast damping of 
LFP autocorrelogram. The lower and the upper layer have now the same average firing rate and 
are on average in an in-phase locking. 

Figure S12. Alternative parameter choices: network with densified inhibition. We assumed 
in the main text that the probability of inhibitory connection is four times larger than the 
probability of excitatory connections. Some Experimental studies like reference [72], however, 
report a probability of inhibitory connection ten times larger than for excitatory connections. We 
show here results for a 1:10 ratio of excitatory to inhibitory connection probability (probabilities 
used are p^^^ = 0.6 and p^^^ = 0.3 for intra-layer inhibitory connections and p^^^ = 0.3 and 
p^^^ = 0.0 for inter-layer inhibitory connections). A: raster plot of the evoked activity of the 
upper layer excitatory population for a 95% level of contrast stimulus. B: autocorrelogram of 
the evoked LFP (for different network sizes). Units are in nA^. Note that synchronous chaos is 
still present, as evidenced by the fast damping of LFP autocorrelogram, accelerating for larger 
network sizes. An additional effect of increased inhibitory density is a stronger tendency to 
resonate for low contrast stimuli. A weak "hump" at frequencies close to 45 Hz is visible even 
in the spectrum of spontaneous activity (not shown). 

Figure S13. E-I connectivity is not required for the generation of oscillations. With our 
parameter choices, oscillations are generated thanks to delayed mutual inhibition. Excitatory 
neurons indeed are not required for the generation of oscillations, but are entrained by the 
oscillation paced by inhibitory cells. This can be proven by numerical simulations in which 
the activity of excitatory neurons is completely suppressed by a strong hyperpolarizing current 
(raster plot in panel A) or in which synapses from excitatory to inhibitory neurons are removed 
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(raster plot in panel B). In both cases the drive to inhibitory neurons in order to maintain their 
rate of activity unchanged. Note that oscillations continue to exist and their frequency does not 
increase consistently. 

Figure S14. Fluctuations for different noise regimes. The fluctuation level of input cur- 
rents can be controlled by acting on the input rates and peak synaptic conductances. Small 
peak coupling conductances and large input rates yield a quasi tonic input ( "small" -variance 
noise). Conversely, stronger peak coupling conductances and smaller input rates give rise to 
input currents with similar average value but stronger fluctuations in time ("large" -variance 
noise). The net input conductance (red = excitatory, blue = inhibitory) of an upper layer ex- 
citatory neuron driven by a full contrast stimulus is shown in panels A (small-variance noise) 
and C (large- variance noise). Subthreshold voltage fluctuation strength is plotted against tuned 
LGN input rate in panels B (small- variance noise) and D (large- variance noise). The rate ranges 
are different for small- and large-variance noises, but are meant to correspond conventionally in 
both cases to the 0%-100% contrast range (see Tables S2 and S3). For both noise regimes, the 
mean excitatory conductance is of ~ 3 — 5nS for the spontaneous activity and of ~ 20nS for 
full contrast stimuli and the mean inhibitory conductance is of ~ 4 — 6nS for the spontaneous 
activity and of ~ 40 — 50nS for full contrast stimuli. At high contrast, however, fluctuations 
in conductance are much stronger for small-variance input noise, because recurrent inputs are 
highly synchronous. Sub-threshold voltage fluctuations at high contrast are comparably strong 
for both noise regimes, because for small-variance noise weaker fluctuations in the input are am- 
plified by strong conductance fluctuations. At low contrast, when the dynamics is asynchronous 
for both noise regimes, voltage fluctuations are stronger for large-variance noise. 

Figure S15. High contrast dynamics for large-variance noise. Dynamics of the upper layer 
for the presentation of a 95%-contrast stimulus. Input noise parameters are reported in Table 
SI. A: raster plot of the excitatory population activity and associated time-histogram of the rate 
of spiking cells. The histogram bar heights denote the fraction of upper layer excitatory cells 
firing in the bin. Bin-size is 2 ms. B: spike trains of 6 excitatory cells highly activated by the 
presented stimulus. C: membrane potential traces for two neurons stimulated simultaneously 
at close-to-preferred orientation (2 top neurons of Panel B in red and green). This dynamics is 
asynchronous, as indicated by the scaling analyses of Figure S16 C-D. 

Figure S16. Temporal decorrelation and spectra of LFPs for large-variance noise. Input 
noise parameters are reported in Table SI. A: autocorrelogram of the LFP evoked by a high 
contrast stimulation. Units are nA^. B: power spectra of evoked LFP for various contrast levels. 
C: scaling with network size of the 95%-contrast synchrony factor x- D: scaling with network 
size of the 95%-contrast LFP autocorrelogram. The dashed line is a power-law with exponent 
-0.5. This scaling is indicative of an asynchronous state. Units are nA^. 

Figure S17. Alternative parameter choices: network with a non modulated spatial profile of 
inter-layer excitation. With the parameter choices assumed in the main text, the integrated effect 
of the inter- layer coupling is inhibitory. It is however moderately excitatory between neurons 
in close vertical alignment, due to the strong spatial modulation of the inter- layer excitation 
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profile. We show here results for the case in which the spatial modulation of inter-layer excitation 
is removed and an equivalent average level of inter-layer excitation is used, but spread across all 
the angular distances (i.e. plnter-iayer E ~ raster plot of the evoked activity of the upper 

layer excitatory population for a 95% level of contrast stimulus. B: autocorrelogram of the 
evoked LFP. Note that synchronous chaos disappears, replaced by almost periodic oscillations, 
very similar to the case of uncoupled layers (F = 0, sec Figure S9). Conversely, removal of 
inter-layer inhibition would further strengthen synchronized chaos (not shown). 

Text SI. Detailed methods for chaos assessment. Section 1: Determination of the minimum 
embedding dimension. Section 2: Extraction of the largest Lyapunov exponent Xmax- 

Text S2. Correspondence between contrast and LGN input rate. Extended description of 
the rationale behind the mapping between contrast and noise input parameters. 

Table SI. Strong noise LGN input parameters. Parameters of the LGN input to the network 
for the high contrast strong noise regime. See Table S3 for more details. 

Table S2. Correspondence between C and Rq^^ for small-variance noise. Correspondences 
are computed approximately, assuming that each cell receives 30 independent AMPA synaptic 
inputs from LGN (see Text S2). For the response of a single LGN cell we assumed ro = 5 Hz 
and r\ = 48 Hz. 

Table S3. Correspondence between C and R^'-^^ for large-variance noise. Correspondences 
are computed approximately, assuming that each cell receives 10 AMPA synapses from 3 inde- 
pendent LGN neurons (see Text S2). For the response of a single LGN cell we assumed ro = 5 
Hz and n = 32 Hz. 
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Table 1. Parameters for model neurons 





Excitatory 


Inhibitory 


' m 


23.30 ms 


11.65 ms 


c 


0.26 nF 


0.13 nF 


vl 


-57.8 mV 


-57.8 mV 




-45.2 mV 


-45.2 mV 


AT 


1.2 ms 


1.2 ms 



Parameters (without time dependency) of model excitatory and inhibitory EIF neurons. 



Table 2. Soft refractoriness parameters 



X 




TA,x 


Bx 


rB,x 


Vl 


22.9 mV 


14.7 ms 


13.5 mV 


76.2 ms 


Vt 


10.0 mV 


17.7 ms 






'\-/Tm 


0.14 ms~^ 


14.3 ms 







Parameters of time-dependent after-spike relaxation of excitatory EIF model neurons 
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Table 3. Synaptic time-constants and efficacies 





Tr (ms) 


Td (ms) 


g (nS) 


PSP (mV) 


AMPA on excitatory 


1 


3 


1.0 


0.84 


AMPA on inhibitory 


1 


3 


1.5 


2.07 


GABA on excitatory 


1 


4 


4.0 


1.13 


GABA on iniiibitory 


1 


2 


4.0 


1.36 


NMDA on excitatory 


3 


80 


0.14 


0.50 



Synaptic parameters for a network of Ne = 4000 neurons and Nj = 1000 neurons: synaptic 
rise (Tr) and decay (r^) times, peak synaptic conductance (g) and peak postsynaptic potential 
PSP. 
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Figure 9 



Table 4. Synaptic latencies 





d (ms) 


Intra-layer synapse 


1.0 


Inter-layer synapse (upper to lower layer) 


3.0 


Inter-layer synapse (lower to upper layer) 


1.0 



Synaptic latencies (d) depending on the relative position of pre- and post-synaptic neurons. 




100 200 
PST (ms) 



i I I 1.1 III I 




CO 



CO 



CO 



100 200 
PST (ms) 



I liL^AA^ 



50 ms 



50 ms 



Figure 10 
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Figure 11 
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Table 5. Probabilities of connection 







pW 


Mean target E cells 


Mean target I cells 


Intra-layer E to E or I 


0.06 


0.06 


~240 


-60 


Intra-layer I to E or I 


0.24 


0.12 


-960 


-240 


Upper E to Lower E or I 


0.06 


0.18 


~240 


~60 


Lower E to Upper E only 


0.06 


0.18 


~240 




Lower E to Upper I only 


0.07 


0.16 




~70 


Inter-layer I to E or I 


0.12 


0.00 


--480 


-120 



Probabilities of connection. The connection probability parameters p^^^ and p^^^ are given for a 
network size of Ne = 4000 and Nj = 1000 per layer. 



Table 6. LGN input 



9LGN 


1 nS 


Ro 


150 Hz 


Ri 


2850 Hz 


ro 


5 Hz 


ri 


48 Hz 



Parameters of the LGN input to the network Rlgn{C). See Text and Tables SIO for more 
details. 



Table 7. Background cortical noise 





10 nS 




10 Hz 


CTbg 


1 Hz 


Ug 


10 ms 



Parameters of the background cortical input R^g. 



